This is an R Markdown Notebook to study boundary layer meteorology.
Create the dataset - - From Great Plains Field Prpgram at O’Neil Nebraska
TatDpoint025 <- c(25.54,24.84,25.77,29.42,33.25,35.25,34.84,32.63,30.07,28.42,27.09,26.09,25.30)
TatDpoint05 <- c(26.00,25.30,25.35,27.36,30.32,32.63,33.20,32.05,30.20,28.74,27.50,26.60,25.83)
TatDpoint10 <- c(26.42,25.84,25.42,25.98,27.62,29.52,30.62,30.62,29.91,28.84,27.84,27.06,26.40)
TatDpoint20 <- c(26.32,25.97,25.56,25.39,25.57,26.11,26.88,27.41,27.68,27.57,27.22,26.87,26.53)
TatDpoint40 <- c(24.50,24.48,24.36,24.34,24.27,24.24,24.26,24.32,24.47,24.64,24.73,24.78,24.84)
datetime <- c('2015-08-31 04:35','2015-08-31 06:35','2015-08-31 08:35','2015-08-31 10:35','2015-08-31 12:35','2015-08-31 14:35','2015-08-31 16:35','2015-08-31 18:35','2015-08-31 20:35','2015-08-31 22:35','2015-09-01 00:35','2015-09-01 02:35','2015-09-01 04:35')
blaygpfp <- data.frame(TatDpoint025,TatDpoint05,TatDpoint10,TatDpoint20,TatDpoint40,datetime)
check the data
#convert the date format
blaygpfp$datetime<- as.POSIXct(blaygpfp$datetime, "%Y-%m-%d %H:%M",tz="US/Pacific")
unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'Warning in as.POSIXlt.POSIXct(x, tz) :
unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'
#TODO Take the annnoying date fotmat error out
str(blaygpfp)
'data.frame': 13 obs. of 6 variables:
$ TatDpoint025: num 25.5 24.8 25.8 29.4 33.2 ...
$ TatDpoint05 : num 26 25.3 25.4 27.4 30.3 ...
$ TatDpoint10 : num 26.4 25.8 25.4 26 27.6 ...
$ TatDpoint20 : num 26.3 26 25.6 25.4 25.6 ...
$ TatDpoint40 : num 24.5 24.5 24.4 24.3 24.3 ...
$ datetime : POSIXct, format:
unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'
"2015-08-31 04:35:00" "2015-08-31 06:35:00" "2015-08-31 08:35:00" ...
Plot the data
plot(blaygpfp)
library(ggplot2)
library(dplyr)
library(reshape)
Attaching package: ‘reshape’
The following object is masked from ‘package:dplyr’:
rename
Plot Temperature waves as function of time and depth.
blaygpfp %>% ggplot() +
geom_point(aes(x=datetime,y=TatDpoint025,colour = as.factor(0.025), shape=as.factor(0))) +
stat_smooth(aes(x = datetime, y = TatDpoint025), se = FALSE) +
geom_point(aes(x=datetime,y=TatDpoint05,colour = as.factor(0.05))) +
stat_smooth(aes(x = datetime, y = TatDpoint05), se = FALSE) +
geom_point(aes(x=datetime,y=TatDpoint10,colour = as.factor(0.10))) +
stat_smooth(aes(x = datetime, y = TatDpoint10), se = FALSE) +
geom_point(aes(x=datetime,y=TatDpoint20,colour = as.factor(0.20))) +
stat_smooth(aes(x = datetime, y = TatDpoint20), se = FALSE) +
geom_point(aes(x=datetime,y=TatDpoint40,colour = as.factor(0.40))) +
stat_smooth(aes(x = datetime, y = TatDpoint40), se = FALSE) +
labs(colour = "Depth") +
xlab("Date") + ylab("Temperature")
# check if line plot is better
With markers
#normalize temperature
blaygpfp %>% ggplot() +
geom_point(aes(x=datetime,y=TatDpoint025, shape=as.factor(0.025))) +
geom_line(aes(x = datetime, y = TatDpoint025)) +
geom_point(aes(x=datetime,y=TatDpoint05,shape = as.factor(0.05))) +
geom_line(aes(x = datetime, y = TatDpoint05)) +
geom_point(aes(x=datetime,y=TatDpoint10,shape = as.factor(0.10))) +
geom_line(aes(x = datetime, y = TatDpoint10),) +
geom_point(aes(x=datetime,y=TatDpoint20,shape = as.factor(0.20))) +
geom_line(aes(x = datetime, y = TatDpoint20)) +
geom_point(aes(x=datetime,y=TatDpoint40,shape = as.factor(0.40))) +
geom_line(aes(x = datetime, y = TatDpoint40)) +
labs(colour = "Depth") +
xlab("Date") + ylab("Temperature")+
ggtitle("Temperatures waves - function of time and depth")
# check if line plot is better
Housekeeping
# reshape it
blaygpfp_melt <- melt(blaygpfp,id="datetime")
Warning in as.POSIXlt.POSIXct(x, tz) :
unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'
# set the depth
blaygpfp_melt$depth <- 0
Warning in as.POSIXlt.POSIXct(x, tz) :
unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'
blaygpfp_melt[blaygpfp_melt$variable=='TatDpoint05',]$depth = 0.05
Warning in as.POSIXlt.POSIXct(x, tz) :
unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'
blaygpfp_melt[blaygpfp_melt$variable=='TatDpoint10',]$depth = 0.10
Warning in as.POSIXlt.POSIXct(x, tz) :
unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'
blaygpfp_melt[blaygpfp_melt$variable=='TatDpoint20',]$depth = 0.20
Warning in as.POSIXlt.POSIXct(x, tz) :
unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'
blaygpfp_melt[blaygpfp_melt$variable=='TatDpoint40',]$depth = 0.40
Warning in as.POSIXlt.POSIXct(x, tz) :
unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'
blaygpfp_melt[blaygpfp_melt$variable=='TatDpoint025',]$depth = 0.025
Warning in as.POSIXlt.POSIXct(x, tz) :
unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'
library(lubridate)
Attaching package: ‘lubridate’
The following object is masked from ‘package:reshape’:
stamp
The following object is masked from ‘package:base’:
date
#In check_tzones(e1, e2) : 'tzone' attributes are inconsistent
blaygpfp_melt$datetime <- ymd_hms(blaygpfp_melt$datetime)
unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'
#Do increment of the time intervals
blaygpfp_melt %>% filter(datetime == ymd_hms("2015-09-01 04:35:00"))
# set the time intervals
blaygpfp_melt$interval <- 0
blaygpfp_melt[datetime == ymd_hms("2015-08-31 06:35:00"), ]$interval<- 2
blaygpfp_melt[datetime == ymd_hms("2015-08-31 08:35:00"), ]$interval<- 4
blaygpfp_melt[datetime == ymd_hms("2015-08-31 10:35:00"), ]$interval<- 6
blaygpfp_melt[datetime == ymd_hms("2015-08-31 12:35:00"), ]$interval<- 8
blaygpfp_melt[datetime == ymd_hms("2015-08-31 14:35:00"), ]$interval<- 10
blaygpfp_melt[datetime == ymd_hms("2015-08-31 16:35:00"), ]$interval<- 12
blaygpfp_melt[datetime == ymd_hms("2015-08-31 18:35:00"), ]$interval<- 14
blaygpfp_melt[datetime == ymd_hms("2015-08-31 20:35:00"), ]$interval<- 16
blaygpfp_melt[datetime == ymd_hms("2015-08-31 22:35:00"), ]$interval<- 18
blaygpfp_melt[datetime == ymd_hms("2015-09-01 00:35:00"), ]$interval<- 20
blaygpfp_melt[datetime == ymd_hms("2015-09-01 02:35:00"), ]$interval<- 22
blaygpfp_melt[datetime == ymd_hms("2015-09-01 04:35:00"), ]$interval<- 24
#'data.frame': 65 obs. of 3 variables:
# $ datetime: POSIXct, format: "2015-08-31 04:35:00" "2015-08-31 06:35:00" #$ variable: Factor w/ 5 levels "TatDpoint025",..: 1 1 1 1 1 1 1 1 1 1 ...
# $ value : num 25.5 24.8 25.8 29.4 33.2 ...
Plot Time waves as function of depth and temperature
#normalize temperature
blaygpfp_melt$scaleT <- scale(blaygpfp_melt$value) %>% as.vector()
# We want for hours 0, 4,8,12,16,20
blaygpfp_melt %>% filter(interval %in% c(0, 4,8,12,16,20)) %>% ggplot() +
geom_point(aes(x=value,y=depth,color=as.factor(interval),shape=as.factor(interval)))+
#geom_line(aes(x=value,y=depth,color=as.factor(interval),shape=as.factor(interval)))+
stat_smooth(aes(x=value,y=depth,color=as.factor(interval),shape=as.factor(interval)))+
xlab("Temperature(deg C)") + ylab("Depth")+
ggtitle("Temperatures at various times - function of z and T")
Ignoring unknown aesthetics: shape
# check if line plot is better
library(ggplot2)
amp <- c(5,3.6,2.5,1.5,.5)
depth <- c(0.025,0.05,.1,.2,.4)
plotampdepthdf <- data.frame(depth,amp)
lm_eqn <- function(df){
m <- lm(amp ~ depth, df);
eq <- substitute(italic(y) == a + b %.% italic(x)*","~~italic(r)^2~"="~r2,
list(a = format(coef(m)[1], digits = 2),
b = format(coef(m)[2], digits = 2),
r2 = format(summary(m)$r.squared, digits = 3)))
as.character(as.expression(eq));
}
p <- ggplot(data = plotampdepthdf, aes(x = depth, y = amp)) +
geom_smooth(method = "lm", se=FALSE, color="black", formula = y ~ x) +
geom_point()
p1 <- p + geom_text(x = 0.2, y = 3, label = lm_eqn(plotampdepthdf), parse = TRUE)
p1
With log Amp
library(ggplot2)
amp <- c(5,3.6,2.5,1.5,.5)
depth <- c(0.025,0.05,.1,.2,.4)
amp <- log(amp)
plotampdepthdf <- data.frame(depth,amp)
lm_eqn <- function(df){
m <- lm(amp ~ depth, df);
eq <- substitute(italic(y) == a + b %.% italic(x)*","~~italic(r)^2~"="~r2,
list(a = format(coef(m)[1], digits = 2),
b = format(coef(m)[2], digits = 2),
r2 = format(summary(m)$r.squared, digits = 3)))
as.character(as.expression(eq));
}
p <- ggplot(data = plotampdepthdf, aes(x = depth, y = amp)) +
geom_smooth(method = "lm", se=FALSE, color="black", formula = y ~ x) +
geom_point()
p1 <- p + geom_text(x = 0.2, y = .5, label = lm_eqn(plotampdepthdf), parse = TRUE)+
xlab("Depth") + ylab("log(amplitude)") + ggtitle("Log of Amplitude")
p1
amp <- c(5,3.6,2.5,1.5,.5)
depth <- c(0.025,0.05,.1,.2,.4)
plotampdepthdf <- data.frame(depth,amp)
poly.model <-lm(amp ~ poly(depth,3,raw=TRUE), df)
summary(poly.model)
Call:
lm(formula = amp ~ poly(depth, 3, raw = TRUE), data = df)
Residuals:
1 2 3 4 5
0.111083 -0.208280 0.121497 -0.026035 0.001736
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.2571 0.4958 12.620 0.0503 .
poly(depth, 3, raw = TRUE)1 -60.9184 13.9633 -4.363 0.1434
poly(depth, 3, raw = TRUE)2 256.3255 89.8197 2.854 0.2146
poly(depth, 3, raw = TRUE)3 -350.0562 146.3810 -2.391 0.2521
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.2668 on 1 degrees of freedom
Multiple R-squared: 0.9943, Adjusted R-squared: 0.977
F-statistic: 57.69 on 3 and 1 DF, p-value: 0.09641
lm_eqn <- function(df){
m <- lm(amp ~ poly(depth,3,raw=TRUE), df);
eq <- substitute(italic(y) == a + b %.% italic(x)*","~~italic(r)^2~"="~r2,
list(a = format(coef(m)[1], digits = 2),
b = format(coef(m)[2], digits = 2),
r2 = format(summary(m)$r.squared, digits = 3)))
as.character(as.expression(eq));
}
p <- ggplot(data = plotampdepthdf, aes(x = depth, y = amp)) +
geom_smooth(method = "lm", se=FALSE, color="black", formula = y ~ x) +
geom_point()
p1 <- p + geom_text(x = 0.2, y = 3, label = lm_eqn(plotampdepthdf), parse = TRUE) +
xlab("Depth") + ylab("Amplitude") + ggtitle("Polynomial fit(Order 3)")
p1
plm_eqn <- function(df){
m <- lm(amp ~ poly(depth,3,raw=TRUE), df);
eq <- substitute(italic(y) == a + b %.% italic(x)*","~~italic(r)^2~"="~r2,
list(a = format(coef(m)[1], digits = 2),
b = format(coef(m)[2], digits = 2),
r2 = format(summary(m)$r.squared, digits = 3)))
as.character(as.expression(eq));
}
predictedamps <- predict(poly.model,list(depth=plotampdepthdf$depth))
p <- ggplot(data = plotampdepthdf, aes(x = depth, y = amp)) +
geom_point()
p1 <- p + geom_line(aes(x=depth,y=predictedamps)) +
geom_text(x = 0.2, y = 3, label = plm_eqn(plotampdepthdf), parse = TRUE) +
xlab("Depth") + ylab("Amplitude") + ggtitle("Predicted with actual")
p1
library(ggplot2)
lag <- c(5,3.6,2.5,1.5,.5)
depth <- c(0.025,0.05,.1,.2,.4)
plotampdepthdf <- data.frame(depth,amp)
lm_eqn <- function(df){
m <- lm(amp ~ depth, df);
eq <- substitute(italic(y) == a + b %.% italic(x)*","~~italic(r)^2~"="~r2,
list(a = format(coef(m)[1], digits = 2),
b = format(coef(m)[2], digits = 2),
r2 = format(summary(m)$r.squared, digits = 3)))
as.character(as.expression(eq));
}
p <- ggplot(data = plotampdepthdf, aes(x = depth, y = amp)) +
geom_smooth(method = "lm", se=FALSE, color="black", formula = y ~ x) +
geom_point()
p1 <- p + geom_text(x = 0.2, y = 3, label = lm_eqn(plotampdepthdf), parse = TRUE)
p1
We need the time at which the maxima occurs, so do the lazy way to get the maxima by installing plotly to get visual cues.
blaygpfp %>% ggplot() +
stat_smooth(aes(x = datetime, y = TatDpoint025,colour = as.factor(0.025)), se = FALSE) +
stat_smooth(aes(x = datetime, y = TatDpoint05,colour = as.factor(0.05)), se = FALSE) +
stat_smooth(aes(x = datetime, y = TatDpoint10,colour = as.factor(0.10)), se = FALSE) +
stat_smooth(aes(x = datetime, y = TatDpoint20,colour = as.factor(0.20)), se = FALSE) +
stat_smooth(aes(x = datetime, y = TatDpoint40,colour = as.factor(0.40)), se = FALSE) +
labs(colour = "Depth(m") +
xlab("Date") + ylab("Temperature (deg C)")+
ggtitle("Soil Temperatures at various depths GPFP (O'Neil Nebraska) ")
# check if line plot is better
ggplotly()
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'
blaygpfp %>% ggplot() +
geom_line(aes(x=datetime,y=TatDpoint025,colour = as.factor(0.025))) +
geom_line(aes(x=datetime,y=TatDpoint05,colour = as.factor(0.05))) +
geom_line(aes(x=datetime,y=TatDpoint10,colour = as.factor(0.10))) +
geom_line(aes(x=datetime,y=TatDpoint20,colour = as.factor(0.20))) +
geom_line(aes(x=datetime,y=TatDpoint40,colour = as.factor(0.40))) +
labs(colour = "Depth(m") +
xlab("Date") + ylab("Temperature (deg C)")+
ggtitle("Soil Temperatures at various depths GPFP (O'Neil Nebraska) ")
# check if line plot is better
ggplotly()
unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'
maxtempat <- c(14,16,16,20,4)
depth <- c(0.025,0.05,.1,.2,.4)
plothourdepthdf <- data.frame(depth,maxtempat)
p <- ggplot(data = plothourdepthdf, aes(x = depth, y = maxtempat)) +
geom_point()
p
Vertical soil temperatures at 0435, 0835, 1235,1635,2035 and 0035
library(lubridate)
fix the timestamp